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ABSTRACT:  Flavins  are  known  to  enhance  extracellular  electron  transfer  (EET) 
in  Shewanella  oneidensis  MR-1  bacteria,  which  reduce  electron  acceptors  through 
outer-membrane  (OM)  cytochromes  c.  Free-shuttle  and  bound-redox  cofactor 
mechanisms  were  proposed  to  explain  this  enhancement,  but  recent  electro¬ 
chemical  reports  favor  a  flavin-bound  model,  proposing  two  one-electron 
reductions  of  flavin,  namely,  oxidized  (Ox)  to  semiquinone  (Sq)  and  semiquinone 
to  hydroquinone  (Hq),  at  anodic  and  cathodic  conditions,  respectively.  In  this 
work,  to  provide  a  mechanistic  understanding  of  riboflavin  (RF)  binding  at  the 
multiheme  OM  cytochrome  OmcA,  we  explored  binding  configurations  at  hemes 

2,  5,  7,  and  10.  Subsequently,  on  the  basis  of  molecular  dynamics  (MD)  simulations,  binding  free  energies  and  redox  potential 
shifts  upon  RF  binding  for  the  Ox/Sq  and  Sq/Hq  reductions  were  analyzed.  Our  results  demonstrated  an  upshift  in  the  Ox/Sq 
and  a  downshift  in  the  Sq/Hq  redox  potentials,  consistent  with  a  bound  RF— OmcA  model.  Furthermore,  binding  free  energy 
MD  simulations  indicated  an  RF  binding  preference  at  heme  7.  MD  simulations  of  the  OmcA— MtrC  complex  interfacing  at 
hemes  5  revealed  a  small  interprotein  redox  potential  difference  with  an  electron  transfer  rate  of  107—  108/s. 


I.  INTRODUCTION 

Microbial  nanowires  have  attracted  great  interest  recently,  for 
example,  for  bioenergy,1  including  for  use  in  microbial  fuel 
cells,"  catalysis/1  reduction  of  graphene  oxide  to  mitigate  use  of 
toxic  chemicals,4  protection  of  steel  from  corrosion,'  or  in 
bioremediation/'  Of  special  interest  is  the  potential  use  of  the 
exoelectrogens  in  neurophysiology  nanoelectronics.7  Indeed,  it 
was  pointed  out  that  exoelectrogens  are  a  critical  part  of 
microbial  electrochemical  technologies  that  aim  to  impact 
applications  in  biosensing  and  biocomputing.1  For  example, 
Leung  et  al.  have  shown  that  Shewanella  oneidensis  MR-1 
exoelectrogens  exhibit  p-type,  tunable  electronic  behavior  with 
mobility  comparable  to  that  of  organic  semiconductors. 
Interestingly,  the  adsorption  and  electron  transfer  on  the  gold 
surface  of  the  deca-heme  cytochrome  MtrF  was  recently 
investigated.  At  the  same  time,  the  extracellular  electron 
transfer  (EET)  for  the  anaerobic  respiration  process  in  S. 
oneidensis  MR-1  that  can  use  multiple  electron  acceptors, 
including  solid  minerals,  was  extensively  investigated, 1 1  where 
reduction  of  electron  acceptors  requires  electron  transfer  across 
the  outer  membrane  (OM).  The  so-called  “electron  conduit”  is 
comprised  of  multiheme  c-type  cytochromes,  and  electrons 
generated  in  the  cell  are  transferred  toward  the  extracellular 
anode.  This  biological  conduit  can  also  transport  electrons  into 
cells  under  cathodic  conditions.1" 

El-Naggar  and  co-workers1'1’14  proposed  that  the  electron 
transport  occurs  through  multistep  hopping  along  chains  of 
hemes  in  MtrCAB— OmcA  complexes.1 1  MtrC  and  OmcA  are 
deca-heme  cytochromes  located  at  the  OM  exterior  surface.17'1' 

'^•ACS  Publications 


MtrB  is  an  OM  /1-barrel  porin  comprising  28  transmembrane  /?- 
strands,  and  MtrA  is  a  deca-heme  periplasmic  protein  which 
associates  with  the  OM  at  the  interior  surface.  '  MtrF,  MtrD, 
and  MtrE  are  homologues  of  MtrC,  MtrA,  and  MtrB. 
MtrC(F)  transfer  electrons  directly  to  an  extracellular  electron 
acceptor  or  indirectly  through  other  OM  cytochromes  such  as 
OmcA.14,1  The  multiheme  OM  cytochromes  are  organized 
into  four  clades,  and  crystal  structures  were  determined  by 
Clarke  and  co-workers  for  MtrF,  UndA,1  OmcA,1’  and 
MtrC.  The  hemes  in  the  OM  cytochromes  are  arranged  as  a 
staggered  cross  with  hemes  2  and  7,  5,  and  10  located  in  the 
opposite  ends  of  shorter  and  longer  crossbeams,  respectively,  as 
shown  in  Figure  1. 

Recently,  Okamoto  et  al.  demonstrated  that  bound  cell- 
secreted  flavin  redox  cofactors  enhance  EET,  1  ~  namely, 
riboflavin  (RF)  and  flavin  mononucleotide  (FMN),  which 
specifically  associate  with  OmcA  and  MtrC,  respectively,  while 
a  free-form  model  of  soluble  flavin  shuttles  was  discussed  by 
Marsili  et  al."4  NMR  spectroscopy  provided  evidence  that 
interactions  between  the  redox  cofactors  with  OmcA  and  MtrC 
occur  near  the  hemes."  This  process  can  provide  ways  to 
enhance  EET,  such  as  in  modulating  the  ionic  strength  to 
strengthen  the  RF  redox  cofactor  binding  affinity," 6  or  through 
development  of  a  synthetic  flavin  biosynthesis  pathway  from 
Bacillus  subtilis  expressed  in  S.  oneidensis.~  El-Naggar  and  co- 
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Figure  1.  Arrangement  of  hemes  in  the  OmcA  OM  cytochrome.11 


workers  also  investigated  the  role  of  the  electrode  type  in 
discerning  EET  of  free-  and  bound-flavin  models  of  S.  oneidensis 
MR-l.2i  Two  distinct  peaks  in  the  differential  pulse 
voltammetry  measurements  were  observed,  corresponding  to 
both  free  and  cytochrome-bound  flavin  models,  due  to  the  use 
of  a  carbon  cloth  electrodes 

An  electrochemical  study1-  has  shown  that,  under  anodic 
conditions,  a  one-electron  reduction  of  the  oxidized  RF  (Ox)  to 
semiquinone  (Sq)  is  involved  in  the  EET  (eq  l),  with  the 
corresponding  redox  potential  (E°)  measured  as  —110  mV 
(against  SHE),  where 

RF(Ox)  +  e“  +  H+  -+  RFH*(Sq)  (l) 

Under  cathodic  conditions,  the  one-electron  reduction  of  RF 
Sq  to  hydroquinone  (Hq)  (eq  2)  has  a  measured  E°  of  —400 
mV,  where 

RFH*(Sq)  +  e"  RFH“(Hq)  (2) 

For  free  RF  in  aqueous  solution,  E°  of  —314  mV  was 
measured  for  the  Ox/Sq  reaction,  while  the  one-electron  Sq/ 
Hq  reduction  has  a  redox  potential  of  —124  mV.-  ’  Therefore, 
for  bound  RF— OmcA  under  anodic  conditions,  the  Ox/Sq 
redox  potential  was  upshifted  by  204  mV.1-  At  the  cell- 
cathode  interface,  E°  of  Sq/Hq  was  downshifted  by  276  mV.  A 
relative  negative  redox  potential  of  the  Sq/Hq  redox  reaction 
for  bound  RF— OmcA  results  in  thermodynamically  unfavorable 
electron  flow  back  from  the  heme  to  the  cathode.12 

However,  despite  these  recent  developments,  details  of  the 
RF -bound  system  were  not  examined,  and  a  crystal  structure  of 
the  complex  is  not  available.  To  provide  an  atomic-level 
understanding  of  the  RF -bound  OmcA  model  and  demonstrate 
consistency  with  experimental  electrochemical  data,  we  report 
an  investigation  of  RF-bound  OmcA  complexes  based  on  free 
energy  molecular  dynamics  (MD)  simulations.  Following 
molecular  docking  of  the  RF  redox  cofactor  at  OmcA  heme 
sites,  which  were  derived  from  the  crystal  structure,1  binding 
free  energies  were  analyzed.  Changes  in  binding  free  energies 
by  mutation  of  hydrophobic  residues  at  binding  sites  were  also 
considered.  On  the  basis  of  the  RF— OmcA  bound  systems, 


redox  potential  shifts  were  calculated  by  MD  simulations  and 
the  predicted  trends  were  found  to  be  consistent  with 
experiment.  Modeling  of  the  complexation  of  OmcA  and 
MtrC  demonstrated  that  the  interprotein  electron  transfer  rate 
is  comparable  to  those  among  hemes  in  MtrF. 

II.  COMPUTATIONAL  METHODS 

11.1.  Docking.  Docking  of  RF  was  carried  out  with 
AutoDock  Vina,’  where  a  semiempirical  free  energy  force 
field  was  used,’1  which  scores  binding  of  a  small  molecule  to 
protein  sites,  including  pairwise  atomic  terms  for  repulsion  and 
London  dispersion,  hydrogen  bonding,  electrostatics,  and 
desolvation.  The  exhaustiveness  parameter  used  was  set  higher 
than  the  default  value  to  ensure  a  better  search.  Polar 
hydrogens  were  added,  and  Kollman’-  charges  were  used  for 
the  protein,  where  the  calculated  electrostatic  potential  is  fit  to 
the  partial  charge  model.  The  input  was  prepared  with  the 
MGL  tools  http://mgltools.scripps.edu/downloads/mgltools- 
l-5-6-release-candidate%20tools.’  The  OmcA  X-ray  structure 
(pdb:  4LMH) 1 1  was  kept  rigid,  but  all  rotatable  bonds  in  RF 
were  included. 

11.2.  MD  Simulations.  Simulations  were  performed  with 
GROMACS  version  5.0.x,34  using  the  AMBER03  force  field  for 
the  protein.  ”  The  general  AMBER  force  field”  (GAFF)  was 
used  for  RF  and  for  dihedral  angle  terms  of  the  hemes.  Force 
constants  for  bonding  and  angular  terms  for  hemes  were 
derived  from  the  Hessian  of  the  optimized  geometry  of  a  model 
compound,’  ,3S  as  shown  in  Figure  SI  in  the  Supporting 
Information.  The  model  system  includes  two  His  ligands 
(modeled  as  methyl-imidazole)  and  heme  B  with  its  covalently 
linked  two  Cys  residues  (CH3CH2S)  through  two  thioether 
bonds  formed  between  the  heme  porphyrin  atoms  CAB  and 
CAC  (labeled  in  Figure  S2a)  with  the  Cys  S  atom.  The 
geometry  of  the  complex  was  optimized  at  the  DFT  M0639/ 6- 
31G*  level  (the  singlet  Fe(ll)  state  for  heme  was  adopted  as 
the  strong  ligand  field,  since  the  bis-imidazole  coordination 
renders  the  ferric  and  ferrous  hemes  to  be  in  low-spin  states16). 
Force  constants  are  listed  in  Table  SI.  RESP32  charges  for  RF 
and  hemes  with  ligated  His  and  bonded  Cys  residues  were 
obtained  by  B3LYP/6-31G*  calculations,  followed  by  analysis 
with  AmberTools  14. 

MD  simulations  at  constant  NPT  for  10  ns  were  then 
performed.  Particle-mesh  Ewald  sums41  were  applied  for 
efficient  treatment  of  long-range  electrostatics  with  a  short- 
range  cutoff  of  1  nm  and  a  relative  strength  of  the  electrostatic 
interaction  at  the  cutoff  of  10-3  at  1  nm  (spacing  of  0.12  nm, 
interpolation  order  of  6).  A  1.2  nm  cutofF  was  used  in 
calculating  the  Lennard-Jones  potential.  The  LINCS  algorithm 
was  chosen  to  reset  bonds  to  their  correct  lengths  after  an 
unconstrained  update.4  The  temperature  was  regulated  by 
velocity  rescaling  with  a  stochastic  term  of  0.1  ps,  and  the 
pressure  by  the  Berendsen  weak  coupling  method.  The  protein 
was  solvated  in  a  box  of  water  with  a  distance  of  1  nm 
separating  the  solute  and  the  box  wall. 

11.3.  Free  Energy  Calculations.  To  calculate  the  free 

energy  change  of  two  different  end  states,  an  alchemical  path 
was  used  to  connect  them.  Stratification  or  multistate  sampling 
was  applied  to  reduce  the  systematic  error,  in  which  a  sequence 
of  intermediate  states  are  introduced  for  the  two  end  states  ( UA 
and  UB),  such  that  U(A)  =  (l  —  A)Ua  +  AUB  and  AG  =  G!  —  G0 
=  !  —  Gj.  To  avoid  numerical  singularity  problems  in 

the  decoupling  procedure,  a  soft-core  potential  was  intro- 
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Figure  2.  Binding  pockets  and  hydrogen  bonding  of  RF  near  hemes  2  (al,  a2),  5  (bl,  b2),  7  (cl,  c2),  and  10  (dl,  d2).  Proximate  residues  and 
hydrogen  bond  distances  (red,  in  A)  are  shown. 


duced,44  where  V£,(r)  =  (l  -  A)VA(RA(r,  X))  +  X VB(RB(r,  X)), 
with  RA(r,  X)  =  {aaAXp  +  r6)1/6  and  RB(r,  X)  =  (atrA6(l  —  X)p  + 
r6)1/6,  and  VA  and  VB  are  van  der  Waals  (vdW)  potentials  for 
states  A  and  B,  respectively,  and  the  soft-core  potential  recovers 
the  A  state  for  X  =  0  and  the  B  state  for  X  =  1.  In  practical 
calculations,  p  =  1,  a  =  0.5,  and  o  =  0.3.  For  each  window  (X,), 
the  following  procedure  was  applied:  (l)  geometry  optimiza¬ 
tion;  (2)  equilibration  at  constant  NVT  for  100  ps;  (3) 
equilibration  at  constant  NPT  for  100  ps;  (4)  MD  simulation 
(sampling)  at  constant  NPT  for  1  ns.  For  the  decoupling 
protocol  in  binding  free  energy  calculations,  5  ns  MD 
simulations  were  carried  out  both  for  NPT  equilibration  and 
for  sampling  to  ensure  sufficient  sampling.  As  X  changes  from  0 
to  1  (see  Table  S2),  the  Bennett  acceptance  ratio  method 
implemented  in  GROMACS  was  used  to  calculate  the  free 
energy  difference  between  different  states.  For  each  window  Xa 
the  potential  energy  difference  in  both  the  forward  A  lfik  = 


Ui+i(xi,k)  ~  and  reverse  Al^.  =  Ui+1(xi+hk)  -  D,(xj+ljt) 

directions  was  sampled.  The  free  energy  difference  A  A  between 
window  i  and  i  +  1  was  computed  by  solving 

\^ni  1 

Zjk=l  1  +  "L  exp[(Auft  -  AA)  /  EbT] 


Y!)  , - w: - 1 - where  n,(nf)  are  the  number  of 

^ k=1  l  +  ^exp[(AA-Auy /kBT]  'V 

sampling  configurations  at  window  i(i  +  l). 


In  binding  free  energy  calculations,  the  virtual  bond 
algorithm  '  was  followed  to  confine  RF  in  the  binding  pocket 
of  OmcA,  using  six  geometrical  parameters,  where  (r0,  d0,  q>0) 
are  the  distance,  angle,  and  dihedral  angle  that  describe  the 
position  of  the  ligand  relative  to  the  protein,  and  (a0,  /?0,  y0)  are 
the  angles  and  dihedral  angles  that  describe  its  orientation 
(defined  in  Figure  S2a).  The  translational  so-called  restraint 
potential  in  polar  coordinates  implemented  in  GROMACS  is 
given  by  u(r,  0,  (p)  =  ^Kr(r  -  r0 )2  +  Ke(l  -  cos (0  -  0O))  + 
^Kq,(<P  ~  cp0f,  and  the  rotational  restraint  potential  by  u(a,  [), 

r)  =  K( 1  -  cos(«  -  «0))  +  \Kp(P  ~  P0)2  +  ^Kr(r  ~  r0)2- 

For  a  given  binding  molecule  (L),  the  configurational  partition 
function  can  be  written  as  ZL  =  Zj*  X  Z”f,  where  Zj*  arises 
from  translational  and  rotational  external  degrees  of  freedom 
(DOFs)  and  ZA'  from  internal  DOFs.  Zf14  =  8^1/  for  a 
nonlinear  molecule.  For  one  molar  standard  state,  V  equals 
1660  A3.  After  applying  the  restraint  potential  to  the  molecule 
(L*),  we  have  ZL *  =  Z“  X  Z“‘,  such  that  Z^?  =  /"  r2  dr  fg  sin 
0d0  fu_„dcpe-u(r’e’<p)/kT  /S  sin  a  da  /^d/9  fn_„dY<Tu{aM,kT. 
Here  T  is  the  temperature  and  k  is  the  Boltzmann  constant. 
Since  there  are  no  cross  interaction  terms  among  the  six 
external  DOFs,  their  contribution  to  the  partition  function  can 
be  calculated  separately  with  numerical  integration  methods 
such  as  Simpson’s  rule." 
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III.  RESULTS  AND  DISCUSSION 

111.1 .  Binding  Free  Energies.  In  MtrF,  the  binding  site  of  a 
soluble  substrate  was  proposed  to  be  near  hemes  2  or  7, 16  in 
UndA  near  heme  7, 18  and  for  OmcA,  the  dimer  structure  was 
reported  to  interface  at  hemes  5,  while  heme  10  was  assumed 
the  electron  egress  site.'  A  molecular  docking  study  of  FMN 
binding  at  MtrC  hemes  found  overall  weak  binding  with  a  slight 
preference  for  heme  2,4"  while  a  similar  investigation  for 
anthraquinone-2, 6-disulfonate  docking  at  UndA  hemes  has 
shown  preference  for  hemes  2  and  7.“  We  considered  RF 
binding  at  hydrophobic  pockets  of  OmcA  near  hemes  2,  5,  7, 
and  10.  Following  identification  of  hydrophobic  pockets  in  the 
OmcA  X-ray  structure,1  molecular  docking  calculations  were 
performed.  Subsequently,  starting  from  the  docked  RF  bound 
OmcA  structures  with  strongest  binding  affinity,  all-atom  MD 
simulations  with  explicit  solvent  were  carried  out.  The  resulting 
RF— OmcA  structures  are  shown  in  Figure  2,  indicating 
hydrogen  bonding  and  hydrophobic  residues  at  the  binding 
sites. 

Examination  of  the  binding  sites  demonstrates  that,  for  RF 
bound  near  hemes  5,  7,  and  10,  the  flavin  isoalloxazine  ring  is 
almost  parallel  to  the  heme  porphyrin  plane,  and  the  nonpolar 
ring  (l,2-dimethylbenzene  ring,  see  RF  structure  in  Figure  S2b) 
approaches  the  heme,  while  the  polar  ring  with  the  — NH  and 
—CO  groups  reaches  out.  For  RF  near  heme  2,  the 
isoalloxazine  ring  is  approximately  perpendicular  to  the  heme 
porphyrin  plane.  RF’s  polar  side  is  close  to  the  heme,  with  the 
— NH  and  one  —CO  groups  forming  H-bonds  with  nearby 
water  molecules,  and  one  —CO  group  forming  an  H-bond  with 
Lysl23.  The  nonpolar  part  of  RF  is  also  exposed  to  the  solvent. 
For  hemes  5,  7,  and  10,  the  flavin  isoalloxazine  ring  is  partially 
inserted  within  the  binding  pocket,  which  involves  His358, 
Tyr374,  and  Tyr517  for  heme  5,  Tyr436,  and  heme’s  porphyrin 
for  heme  7,  and  Tyr688,  Tyr622,  and  heme  9’s  porphyrin  for 
heme  10.  For  heme  2,  the  flavin  ring  stacks  against  the  Tyr301 
phenol  ring.  In  all  cases,  one  or  more  Tyr  residues  are  involved 
in  the  binding  pocket,  as  is  known  for  flavoproteins,  where 
aromatic  stacking  between  the  flavin  isoalloxazine  ring  and 
hydrophobic  residues,  such  as  Tyr  and  Trp,  plays  a  major  role 
in  the  binding  process.4’  The  ribityl  —OH  groups  (2'-OH,  3'- 
OH,  4'-OH,  and  5-OH  in  Figure  S2b)  form  hydrogen  bonds 
with  the  protein  or  solvent  water  molecules.  For  heme  2,  2'- 
OH,  3'-OH,  and  4'-OH  form  hydrogen  bonds  with  a  nearby 
water  molecule,  while  5-OH  with  heme  2’s  porphyrin  acetate 
group.  For  heme  5,  4'-OH  and  5-OH  form  hydrogen  bonds 
with  a  nearby  water  molecule,  while  2'-OH  and  3-OH  with 
heme  2’s  porphyrin  acetate  group.  For  heme  7,  2'-OH  forms  a 
hydrogen  bond  with  the  Tyr577  —OH  group  and  the  backbone 
—CO  group,  3-OH  with  the  N1  atom,  and  4'-OH  with 
Ser440’s  —OH  group  and  the  backbone  —CO  group  of  Tyr460, 
and  5-OH  with  a  nearby  water  molecule.  For  heme  10,  2'-OH 
forms  a  hydrogen  bond  with  —OH  of  Tyr688,  3  -OH  with  a 
nearby  water  molecule,  and  4  -OH  with  Lys632  and  S'-OH 
with  the  Glu614  acetate  group. 

Binding  free  energies  of  oxidized  RF  near  hemes  2,  5,  7,  and 
10  in  OmcA  were  calculated  by  the  alchemical  double 
decoupling  method  ' 44  with  a  restraint  potential,  as  illustrated 
in  Figure  3.  The  binding  free  energy  (AG^)  of  RF  to  OmcA  in 

agJ 

aqueous  solution  is  described  by  (OmcA)H20  +  (RF)H20 - > 

(OmcA  -  RF)Hi0.46  Following  common  practice,  the 
interaction  between  RF  and  the  environment  for  both  free 


Figure  3.  Thermodynamic  cycle  used  to  calculate  the  RF  binding  free 
energy  at  OmcA  hemes.  AG1  denotes  the  free  energy  change  upon 
turning  off  the  interaction  between  free  RF  and  water;  AGlla  is  the  free 
energy  change  after  exerting  a  restraint  potential  to  confine  the 
position  and  orientation  of  bound  RF;  AGIIb  is  the  free  energy  change 
when  decoupling  the  interaction  between  restrained  bound  RF  (RF*) 
and  the  environment;  AGllc  denotes  the  free  energy  change  of 
removing  the  external  restraint  potential. 


RF  in  solution  and  RF  bound  at  OmcA  were  decoupled,  such 

✓  x  AG1  .  .  .  AGn 

that  (RF)H20  — >  (RF)vac  and  (OmcA  -  RF)H20  - > 

(OmcA)H20  +  (RF)^  and  AGjJ  =  AG1  —  AG11.  A  restraint 
potential  was  applied  to  confine  bound  RF’s  position  and 
orientation  and  avoid  its  wandering  when  interaction  with  the 

AGIIa 

environment  was  turned  off,4'  so  that  (OmcA  —  RF)H20 - > 

AGIIb  agUc 

(OmcA  -  RF*)H20  - »  (OmcA)Hi0  +  (RF*)vac  - > 

(OmcA)H20  +  (RF)vac.  AGjJ  =  AG1  -  AGIIa  -  AGIIb  -  AGUc, 
where  denotes  a  restrained  state.  AG1,  AGUa,  and  AGIIb  were 
readily  calculated  by  free  energy  MD  simulations  using  the 
alchemical  path  and  the  free  energy  change  of  switching  off  the 

(see  Computa- 

/ 

tional  Methods). 

Calculated  binding  free  energies  of  RF  bound  at  wild-type 
(WT)  and  site-mutated  OmcA  with  all  heme  groups  reduced 
are  listed  in  Table  1.  RF  bound  near  hemes  5  and  7  has  a 
relatively  stronger  affinity,  with  binding  free  energies  of  —35.7 
and  —38.2  kj/mol,  respectively,  while  the  binding  affinity  near 
hemes  2  and  10  is  relatively  weaker,  with  binding  free  energies 
of  —20.7  and  —22.6  kj/mol,  respectively,  smaller  by  about  15 
kj/mol.  Interestingly,  the  binding  affinity  of  RF  near  heme  7 
was  —20.5  kj/mol  when  all  heme  groups  were  oxidized, 
decreasing  by  about  18  kj/mol  compared  to  the  respective 
system  with  reduced  hemes.  Note  that  bound  redox  flavin 
cofactors  were  only  observed  for  reduced  MtrC."  ’ 

In  examining  contributions  of  electrostatic  and  vdW 
interactions  to  the  binding  affinity,  we  note  that  the  relatively 
weaker  binding  at  hemes  2  and  10  is  mainly  due  to  a  smaller 
contribution  of  vdW  interactions  because  the  isoalloxazine  ring 
of  RF  at  heme  2  is  exposed  to  the  solvent,  leaving  only  Tyr301 
for  aromatic  stacking.  For  heme  10,  the  RF  isoalloxazine  ring  is 
loosely  stacked  with  neighboring  hydrophobic  residues,  while 
weaker  binding  of  RF  at  heme  7  for  oxidized  OmcA  is  due  to  a 
weaker  electrostatic  interaction.  In  exploring  details  of  binding 


restraint  potential  by  AG  c  =  —kT  In 
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Table  1.  Binding  Free  Energies  (kj/mol)  with  OmcA  Flemes  in  the  Reduced  State'1 


-AGIIb 


binding  site 

AG1 

<1 

1 

electrostatic 

vdW 

total 

<1 

1 

AGg 

2(WT) 

106.04 

-19.31 

-110.88 

-35.61 

-146.49 

39.10 

-20.66 

5(WT) 

106.04 

-12.12 

-95.60 

-74.45 

-170.05 

40.45 

-35.68 

7(WT) 

106.04 

-8.66 

-106.74 

-68.52 

-175.26 

39.69 

-38.19 

7(WT)b 

106.04 

-14.97 

-79.9 

-72.15 

-152.05 

40.49 

-20.49 

7(A436) 

106.04 

-8.59 

-99.05 

-53.29 

-152.34 

39.12 

-15.77 

7(W436) 

106.04 

-10.94 

-94.82 

-65.01 

-159.83 

39.98 

-24.75 

10(WT) 

106.04 

-8.47 

-108.29 

-51.37 

-159.66 

39.46 

-22.63 

10(A688) 

106.04 

-7.77 

-112.31 

-50.85 

-163.16 

39.97 

-24.92 

10(W688) 

106.04 

-12.20 

-107.77 

-57.63 

-165.40 

39.50 

-32.06 

aA  description  of  the  free  energy  terms  is  given  in  the  text.  ^Hemes  in  oxidized  state. 

upon  mutation,  we  found  that  Tyr436/Ala436  and  Tyr436/ 
Trp436  mutations  within  the  binding  pocket  near  heme  7 
decreased  the  binding  affinity  by  22.4  and  13.4  kj/mol, 
respectively,  as  expected.  On  the  other  hand,  Trp688/Ala688 
and  Tyr688/Trp688  mutations  increased  the  binding  affinity 
for  the  binding  pocket  near  heme  10  by  ca.  2.3  and  9.4  kj/mol, 
respectively.  To  assess  the  sensitivity  of  the  binding  affinity  to 
the  details  of  the  binding  site,  alternative  binding  sites  at  hemes 
2,  5,  7,  and  10  were  explored  (see  description  and  Figure  S3  in 
the  Supporting  Information),  resulting  in  binding  energies 
comparable  to  those  reported  here. 

111.2.  Redox  Potential  Shifts.  Next,  we  calculated  the  redox 
potential  shift  of  RF  bound  at  OmcA  compared  to  free  RF  in 
solution.  First,  on  the  basis  of  the  thermodynamic  cycle  shown 
in  Figure  4,  redox  reaction  1  was  considered.  The  difference  in 


Figure  4.  Thermodynamic  cycle  used  to  calculate  the  redox  potential 
shift  when  moving  the  RF  redox  cofactor  from  water  to  the  OmcA 
protein  for  the  one-electron  reduction  reaction  Ox/Sq. 


the  Ox/Sq  reaction  free  energy  can  be  derived  from  the  free 
energy  difference  between  Ox(RF)  and  Sq(RFH)  in  protein 
and  in  water,  where 

AAG  =  AGW^P(RF  +  e“  +  H+  ->  RFH) 

=  AGp,(RF  -  RFH)  -  AG*,(RF  ->  RFH)  (3) 

To  calculate  the  solvation  free  energy  difference  of  RFH  and 
RF  (RFH  has  an  extra  hydrogen  bonded  to  the  N5  atom,  as 
compared  to  RF;  see  the  structure  of  RF  in  Figure  S2b),  we 
designed  reference  states  RFHO  and  RFO,  which  have  the  same 


molecular  structure  as  RFH  and  RF,  respectively,  but  no 
electrostatic  interaction  with  the  environment,  so  that 

AGspjw)(RF  -+  RFH) 

=  AGp(!W)(RF  -»•  RFO)  +  AGp(!W)(RF0  ->  RFHO) 

+  AGp(!W)(RFH0  RFH) 

=  AGp(!W)(RF  -»•  RFO)  -  AGp(!W)(RFH  ->  RFHO) 

-  AGsP|W)(RFH0  -+  RFO)  (4) 

Specifically,  to  change  RF  to  RFO  or  RFH  to  RFHO,  the 
electrostatic  interaction  between  RF  and  the  environment  is 
decoupled,  while,  from  RFHO  to  RFO,  the  vdW  interaction 
between  H(N5)  and  the  environment  is  decoupled,  which  is 
negligible  in  our  calculations. 

Next,  we  considered  reaction  2.  The  redox  potential  shift 
when  moving  RF  from  solution  to  OmcA  for  the  Sq/Hq  redox 
reaction  is  related  to  the  free  energy  difference 


AAG  =  AGP(RFH  +  e-  -+  RFH-) 

-  AGW(RFH  +  e-  ->  RFH-)  (5) 


In  practice,  an  additional  free  RFH  was  placed  and  restrained 
around  2.5  nm  away  from  the  protein  surface.  The  free  energy 
change  upon  electron  transfer  from  free  RFH-  in  the  solvent  to 
RFH  bound  at  OmcA  was  calculated  using  an  alchemical  path, 
in  which  state  A  corresponds  to  RFH-(free)  +  RFH(bound) 
and  state  B  corresponds  to  RFH(free)  +  RFH- (bound).  The 
redox  potential  shift  for  eqs  3  and  5  is  given  by 


AAG 

F 


-AAG  X  10.4 


(6) 


where  the  redox  potential  is  in  mV  and  AAG  in  kj/mol.  Since 
the  redox  potentials  for  the  Ox/Sq  and  Sq/Hq  reactions  were 
measured  under  anodic  and  cathodic  conditions,  oxidized  and 
reduced  heme  states  were  assumed,  respectively. 

The  calculated  results  are  summarized  in  Table  2.  The 
solvation  free  energies  for  RF  and  RFH  are  similar  in  water  yet 
demonstrate  significant  changes  in  the  bound  state,  attributed 
to  differences  in  atomic  partial  charge  distributions  within  the 
protein  environment.  An  upshift  of  the  Ox/Sq  reduction 
potential  of  bound  RF  vs  free  RF  was  indeed  reproduced, 
although  the  calculated  values  of  58,  58,  101,  and  59  mV  for 
binding  sites  near  hemes  2,  5,  7,  and  10  are  smaller  than  the 
experimentally  derived  estimate  of  204  mV.  The  downshifts  of 
the  one-electron  Sq/Hq  for  binding  sites  near  hemes  2,  5,  7, 
and  10  were  calculated  as  —105,  —240,  —469,  and  —498  mV, 
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Table  2.  Free  Energy  Differences  (kj/mol)  for  Reactions  1  and  2“ 


heme  2 

heme  5 

heme  7 

heme  10 

Exp.12 

AGpol(RF  -  RFH) 

Ox/Sq  (eq  l) 

-5.3 

-5.3 

-9.4 

-5.4 

AG",(RF  -  RFH) 

0.3 

0.3 

0.3 

0.3 

AGspol(RF  -  RFH)  -  AG“,(RF  RFH) 

-5.6 

-5.6 

-9.7 

-5.7 

£°  -  e° 

cp  CW 

58 

58 

101 

59 

204 

AGP(RFH  +  e“  ->  RFH-1)  -  AGW(RFH  +  e“  ->  RFH-1) 

Sq/Hq  (eq  2) 

10.2 

23.2 

45.3 

48.0 

0  0 

F  —  F 
cp  cw 

-105 

-240 

-469 

-497 

-276 

is  the  shift  in  reduction  potentials  for  the  Ox/Sq  and  Sq/Hq  redox  reactions  when  RF  changes  from  free-form  in  solution  to  hound  at 
OmcA  (mV). 


respectively,  while  the  experimental  value  is  —276  mV. 
Encouragingly,  the  trend  in  redox  potentials  is  predicted 
correctly,  while  a  quantitative  comparison  to  experiment  for 
this  complex  electrochemical  system  is  beyond  the  scope  of 
current  methodologies  and  available  experimental  data.  For 
example,  the  measured  potential  of  the  two-electron  free-form 
molecular  flavin  reduction  of  —240  mV* 5  is  larger  than  the 
value  of  —207  mV  measured  without  the  cellular  context. 

111.3.  Electron  Transfer  between  MtrC  and  OmcA.  To 
further  understand  EET  in  this  case,  complexation  of  OmcA 
and  MtrC  with  hemes  5  at  the  interface  was  investigated. 
Edwards  et  al.1  pointed  out  that  heme  5  shows  the  largest 
variation  among  crystal  structures,  suggesting  that  this  could  be 
the  site  of  electron  exchange.  In  addition,  it  was  demonstrated 
that  the  OmcA  crystal  structure  and  oligomeric  structure  in 
solution  are  in  excellent  agreement.  OmcA  and  MtrC" '  were 
first  docked  with  ZDOCK,  followed  by  a  MD  simulation  for 
100  ns  at  constant  NPT  (T  =  300  K),  in  which  only  heme  5  of 
OmcA  is  in  the  reduced  state  and  all  other  19  hemes  are  in  the 
oxidized  state.  The  resulting  complex  is  shown  in  Figure  5, 


Figure  S.  Complex  of  OmcA  and  MtrC.  The  distance  between  hemes 
5  is  given  in  A. 


revealing  shape  compensation  at  the  interface  between  OmcA 
and  MtrC  upon  binding,  and  a  close  contact  in  the  region 
around  hemes  7  and  heme  5  with  an  edge-to-edge  distance 
between  the  porphyrins  of  hemes  5  of  5.97  A.  This  is  consistent 
with  a  value  of  5.97  A  calculated  for  the  MtrF— OmcA 
complex,'"  but  our  MtrC— OmcA  complex  exhibits  a  more 


compact  structure.  The  rate  of  electron  transfer  (fcET)  based  on 
Marcus  theory53  is  given  by  log(kET)  =  15  —  0.434/1R  — 
3.1  [(AG  +  X)2/X\,  where  AG  and  A  are  the  Gibbs  free  energy 
and  reorganization  energy  (eV),  respectively.  A  ft  value  of  1.4 
A-1  is  often  assumed,54  intermediate  between  doped  semi¬ 
conductors  and  a  vacuum,55  and  R  is  the  edge— edge  distance 
between  two  heavy  atoms  in  hemes. 

Free  energy  MD  simulations  to  examine  electron  transfer 
from  OmcA’s  heme  5  to  MtrC’s  heme  5  were  performed  using 
an  alchemical  path,  where  states  A  and  B  designate 

HemeS, OmcA*  Fe(Hl)Heme5,MtrC  an<3  Fe(Hl)Heme5,OmcA/  Fe" 

(H)Heme5,MtrO  respectively.  The  reorganization  energy  was 
calculated  using  the  linear  response  approximation  by  A  = 
0-5  ((VB  -  VA)A  -  (VB  -  VA)B),  where  VA  (VB)  are  potential 
energies  when  the  active  site  is  in  state  A  (B),  calculated  by  MD 
simulations.  The  sampling  procedure  is  similar  to  the  reaction 
free  energy  calculation  except  that  only  the  initial  and  final 
states  are  involved  without  intermediate  states.  A  small  reaction 
free  energy  of  about  99  mV  and  reorganization  energy  of  1.1  eV 
were  calculated.  Calculated  uphill  and  downhill  electron 
transfer  rates  (fcEX)  were  2.3  X  107  and  3.9  X  108/s, 
respectively,  comparable  to  the  electron  transfer  rate  between 
intrahemes  in  MtrF.  Note  that  in  vivo  OmcA  forms  a  complex 
with  MtrC  at  the  end  of  the  MtrABC  complex  and  not  with 
free  MtrC. 

IV.  CONCLUSION 

In  this  work,  we  investigated  binding  affinity  and  changes  in  the 
redox  properties  of  riboflavin  upon  binding  to  OmcA  near 
hemes  2,  5,  7,  and  10.  Molecular  docking  combined  with  MD 
simulation  refinement  and  subsequent  MD  free  energy 
calculations  were  performed,  where  derivation  of  force  field 
parameters  for  the  hemes  was  carried  out.  The  calculated 
reduction  potential  was  upshifted  for  the  Ox/ Sq  redox  reaction 
as  compared  to  free  RF  in  solution,  favoring  withdrawal  of 
electrons  from  the  cell,  and  downshifted  for  the  Sq/Hq  redox 
reaction,  consistent  with  experimental  measurements.  Binding 
free  energies  demonstrated  a  binding  preference  of  RF  at 
hemes  5  and  7.  However,  heme  5  may  interface  with  other 
proteins,  thus  leaving  heme  7  as  a  likely  binding  site,  where  the 
binding  affinity  was  also  sensitive  to  Y436/A436  mutation. 
Although  heme  7  is  relatively  insulated  inside  the  protein,  the 
bound  riboflavin  cofactor  extends  toward  the  surface.  This 
observation  could  prompt  experimental  investigation  of 
mutated  OmcA  to  assess  the  role  of  heme  7.  MD  simulations 
for  the  OmcA  and  MtrC  complex  with  interfacing  hemes  5 
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resulted  in  a  relatively  small  edge-to-edge  distance  and  a  small 
redox  potential  difference,  with  kET  comparable  to  that  between 
intraprotein  hemes. 
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